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Abstract 


Despite their promise and ubiquity, Gaussian processes (GPs) can be difficult to use 
in practice due to the computational impediments of fitting and sampling from them. 
Here we discuss a short R package for efficient multivariate normal functions which uses 
the Repp and ReppEigen packages at its core. GPs have properties that allow standard 
functions to be sped up; as an example we include functionality for Toeplitz matrices whose 
inverse can be computed in 0{n^) time with methods due to Trench and Durbin (Golub 
& Van Loan 1996), which is particularly apt when time points (or spatial locations) of a 
Gaussian process are evenly spaced, since the associated covariance matrix is Toeplitz in 
this case. Additionally, we include functionality to sample from a latent variable Gaussian 
process model with elliptical slice sampling (Murray, Adams, & MacKay 2010). 
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1. Introduction 


Many methodologies involving a Gaussian process rely heavily on computationally expensive 
functions such as matrix inversion and the Gholesky decomposition. Rather than create a 
package to solve a particular high-level Gaussian process (GP) task (e.g., expectation prop¬ 
agation, variational inference, regression and classification (Neal 1998)), the aim of FastGP 
is to improve the performance of these fundamental functions in order to help all researchers 
working with GPs. While there exist R packages for sampling from a multivariate normal 
distribution (MVN) or evaluating the density of an MVN, notably MASS and mvtnorm on 
GRAN (Genz & et al. 2014; Venables et al. 2002), we have found such packages can be slow 
in the context of GPs, partially due to unnecessary checks for symmetry and positive definite¬ 
ness (which hold for GPs with commonly used kernels such as squared exponential or Matern 
(Rasmussen & Williams 2005)) or not accounting for the structure (e.g. Toeplitz) of the un¬ 
derlying covariance matrix. Hence, we write functions optimized with Repp and ReppEigen 
(Bates & Eddelbuettel 2013; Eddelbuettel 2013) to make these tasks more computationally 
efficient, and demonstrate their efficiency by benchmarking them against built-in R functions 
and methods from the MASS and mvtnorm libraries. Additionally, we include functionality 
to sample from the posterior of a Bayesian model for which the prior distribution is multi- 
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variate normal using elliptical slice sampling, a task which is often used alongside GPs and 
due to its iterative nature, benefits from a C++ version (Murray, Adams, & MacKay 2010). 


To elaborate, a Gaussian process (GP) is a collection of random variables (i.e., a stochas¬ 
tic process) (Xt) such that any finite subset of these random variables has a joint multivariate 
normal distribution (Grimmet & Stirzaker 2001). Such processes have been used extensively 
in recent decades, particularly in machine learning, spatial and temporal statistics, and com¬ 
puter experiments (Rasmussen & Williams 2005). However, GPs can be difficult to work 
with in practice because they are computationally onerous; to be precise, the density of a 
multivariate normal (MVN) vector in M"' is: 


f(x) 
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This involves the computation of a determinant and inverse of a matrix in generally 

taking O(n^) operations to complete, and the computation of the Cholesky decomposition is 
typically a prerequisite for sampling from a multivariate normal, which also takes O(n^) time 
to complete (Golub X Van Loan 1996). Many spatial models, including those which tackle 
nonstationarity as in Bornn et al. (2012), parametrize the covariance matrix S, and hence for 
Monte Carlo-based inference require repeated recalculations of and |S|. 


2. Key functions and benchmark results 


2.1. Key functions and package organization 

The core functions of FastGP can be categorized into three sets. The first set of functions 
are matrix operations that are necessary for sampling from and evaluating the density of a 
multivariate normal random variable. These are: a function for inversion using RcppEigen, 
a function for inverting a symmetric positive definite Toeplitz matrix in O(n^) time (which, 
as aforementioned, can be useful for inverting a covariance matrix in which the underlying 
points are evenly spaced (Storkey, A.J. 1999)) which uses methods due to Trench and Durbin 
(Golub & Van Loan 1996) written in Repp, and a function for evaluating the Cholesky 
decomposition of a matrix using RcppEigen. To be as explicit as possible, the inversion and 
Cholesky decomposition come directly from RcppEigen (Bates & Eddelbuettel 2013). The 
second set of functions are those which directly simulate from and evaluate the log density of 
a multivariate normal, both in the general case and when the underlying covariance matrix 
is Toeplitz. The final major function included in the package is the elliptical slice sampling 
algorithm for simulating from the posterior of a Bayesian model in which the prior is jointly 
multivariate normal, which tends to outperform classical methods such as Metropolis-Hastings 
computationally, as is evidenced empirically by Murray, Adams, & MacKay (2010). 

2.2. Benchmark results 

Here we use the rbenchmark (Eugster X Leisch 2008) package to demonstrate the efficacy of 
these methods. In particular we test these functions with a mock covariance matrix from a 
square exponential Gaussian process on 200 evenly spaced time points with a and (p arbitrarily 
set to 1 (and hence the covariance matrix is Toeplitz in this case). 
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Table 1: Benchmarking the runtime for functions included in FastGP. The numbers indicate 
how many times faster the functions performed using our package versus standard R, MASS, 
and mvtnorm functions, using the benchmark function from the rbenchmark package. 


FastGP Function 

Standard Function 

Relative Speed Improvement 

rcppeigen_invert_matrix 

solve 

X3.287 

tinv 

solve 

X14.374 

rcppeigen_get_det 

det 

xl.851 

rcpp_log_dmvnorm, istoep= TRUE 

dmvnorm 

xl.966 

rcpp_log_dmvnorm, istoep= FALSE 

dmvnorm 

X.6592 

rcpp_rmvnorm 

rmvnorm 

X23.462 

rcpp_rmvnorm 

mvrnorin 

X22.812 


2.3. Demonstration of elliptical slice sampling 

We consider a model where we observe a “warped” signal s = Asin((t + w)/T) + e where w 
is drawn according to a 0 mean GP with squared exponential kernel, e is drawn according to 
a normal distribution with 0 mean and a = .001, and A and T are known constants. Our 
objective is to perform inference on the latent warping w, and we can do this with the elliptical 
slice sampling function included with FastGP, as in FastGPdemo. r. The function implements 
the algorithm as described by Murray, Adams, & MacKay (2010) and benefits from the use 
of the optimized functions contained within FastGP since each iteration requires several log- 
likelihood evaluations (which may require the evaluation of the log-density of a multivariate 
normal distribution) and drawing from a multivariate normal distribution. This results in 
the following posterior draws for w illustrated in Figure 1. Additionally we benchmark an 
elliptical slice sampler with FastGP functions versus an elliptical slice sampler with standard 
functions below. 

Table 2: Benchmarking the runtime for elliptical slice sampling using functions from FastGP 
versus elliptical slice sampling using standard functions from R, mvtnorm, and MASS. 


FastGP Function 

Standard Function 

Relative Speed Improvement 

rcpp_ess 

standard_ess 

xl.508 


3. Conclusion 

To summarize, we have written an R package FastGP using Repp and ReppEigen for handling 
multivariate normal distributions in the context of Gaussian processes efficiently. Additionally 
we have included functionality to perform Bayesian inference for latent variable Gaussian 
process models with elliptical slice sampling (Murray, Adams, & MacKay 2010). 
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Figure 1: Top: The red line indicates ground truth, the black line indicates the mean of the 
posterior latent samples, and the dotted lines indicate the 2 standard deviations above and 
below the mean, respectively, of the latent samples. Remaining rows show the 100th, 300th, 
500th, 700th, 900th, and 1000th MCMC samples for the warping w respectively, in black, 
compared to ground truth in red. 
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